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ABSTRACT 

The  following  model  problem  examines  the  thermocapillary  feedback  mechanism  im- 
portant at  the  edge  of  weld  pools  and  other  materials  processes.  A  pool  of  liquid  with  a  flat 
horizontal  free  surface  is  bounded  on  one  side  by  a  vertical  solid  wall,  which  is  maintained 
at  a  cold  temperature  to  unit  depth,  and  at  a  warmer  temperature  below;  far  away  the 
fluid  is  at  the  warmer  temperature.  Surface  tension  is  a  decreasing  function  of  tempera- 
ture, so  that  the  surface  thermal  gradient  drives  flow  toward  the  corner.  When  convection 
is  vigorous,  the  flow  compresses  the  thermal  gradient  which  is  driving  the  flow:  this  posi- 
tive feedback  results  in  small  local  length  scales  and  high  velocities  near  the  corner.  This 
problem  is  examined  through  a  detailed  scaling  analysis  and  through  numerical  simulation 
for  a  range  of  parameters. 


THERMOCAPILLARY  FLOW  NEAR  A  COLD  WALL 

1.  INTRODUCTION 

In  the  processing  of  materials,  often  material  is  melted  and  resolidified.  Several  practi- 
cal processes,  e.g.,  welding,  float-zone  purification,  and  Czochralski  crystal  growth,  involve 
a  pool  of  molten  metal  with  a  free  surface,  with  strong  temperature  gradients  along  the 
surface.  Convection  in  the  molten  metal  is  typically  vigorous  and  significant  to  the  results 
of  the  process,  in  that  it  affects  the  size  and  shape  of  the  pool,  the  heat  transfer,  the 
mixing  of  solutes,  and  ultimately  the  microstructure  of  the  finished  product.  The  forces 
driving  the  convection  include  the  variation  of  surface  tension  with  temperature  along 
the  surface  (thermocapillary  forces),  buoyancy  forces  due  to  thermal  (and/or  solutal)  ex- 
pansion, and  electromagnetic  forces  in  the  case  of  arc  welding  or  electron  beam  welding. 
However,  in  many  cases  (e.g.  laser  welding)  thermocapillary  forces  predominate,  and  even 
in  cases  where  other  forces  are  stronger  overall,  there  are  still  important  regions  where  the 
thermocapillary  forces  are  dominant  (i.e.,  cold  corner  regions:  see  Chen.  1987). 

Consequently,  there  have  been  many  theoretical  studies  of  thermocapillary  flows,  pri- 
marily numerical,  and  a  few  analytical  (reviews  are  given  by  Ostrach.  1982,  and  Davis, 
1987).  Cowley  and  Davis  (1983)  analyzed  the  (two-dimensional)  thermocapillary  flow  near 
a  hot  wall  for  vigorous  flow  (large  Marangoni  number):  here  the  fluid  flows  up  the  wall 
then  turns  and  flows  away  along  the  free  surface,  so  this  would  be  called  the  hot  corner 
problem.  The  numerical  studies  of  Zebib  et.  al.  (1985)  of  flow  in  a  rectangular  pool  (2-D) 
with  one  hot  and  one  cold  wall,  however,  show  that  for  moderate  to  small  Prandtl  numbers 
(e.g..  metals)  the  cold  corner  region  has  by  far  the  strongest  effect  on  both  the  flow  and  the 
heat  transfer.  This  result  gives  a  different  overall  scaling  than  that  of  Cowley  and  Davis. 
although  their  local  hot-corner  scaling  was  validated.  Other  numerical  studies  (e.g.,  Zehr 
et.  al..  1987).  when  a  sufficiently  fine  mesh  is  used,  show  similar  strong  flow  at  the  cold 
corner,  where  the  flow  along  the  free  surface  toward  the  cold  wall  compresses  the  thermal 
gradient,  thereby  enhancing  the  flow  and  the  heat  transfer.  Great  care  is  necessary  to 
insure  that  the  small  length  scales  of  this  corner  region  are  resolved  numerically;  this  is 
not  always  the  case  (as  noted  by  Chen,  1987). 

Therefore,  it  is  imperative  to  develop  a  theoretical  understanding  of  the  dynamics 
of  the  cold  corner  region.  Being  a  region  of  intense  heat  transfer,  the  details  of  the  flow 
can  affect  the  shape  of  the  melt  pool  and  the  cooling  rate,  thus  the  microstructure,  of 
the  material.  At  the  least,  the  dependence  of  the  length,  velocity  and  thermal  scales 
on  the  parameters  (Marangoni  number.  Prandtl  number,  Capillary  number)  needs  to  be 
understood  in  order  for  realistic  numerical  models  to  be  designed  in  a  way  to  resolve  the 
details  in  this  important  region.  But  as  yet,  such  understanding  is  lacking.  In  fact,  in  a 
recent  review,  M.  M.  Chen  (1987,  p. 552)  states,  "It  would  seem  then  that  the  structure  of 
the  cold  corner  flow  is  one  of  the  most  critical  issues  to  be  studied  in  the  future." 

To  analyze  the  behavior  of  the  cold  corner  region  without  all  the  complications  of  the 
complex  geometry,  phase  change,  and  time  dependence  inherent  in  real  materials  process- 
ing applications,  a  simplified  model  problem  will  be  considered,  much  like  that  of  Cowley 
and  Davis  (1983),  as  follows.  A  pool  of  liquid  has  a  horizontal  free  surface  ending  at  a 
vertical  wall,  and  the  upper  section  of  the  wall  is  cooled;  the  resulting  thermal  gradient 
drives  thermocapillary  flow  towards  the  cold  corner.  The  depth  and  width  of  the  pool  are 


assumed  large  compared  to  all  local  length  scales  (which  is  reasonable  for  practical  situa- 
tions with  high  Marangoni  numbers),  so  the  pool  appears  semi-infinite  both  horizontally 
and  vertically. 

This  simplified  problem  is  still  complicated,  and  contains  most  of  the  features  of  the 
cold  corner  regions  in  practical  processes,  e.g.,  welds.  The  missing  features  are  phase 
change  and  surface  deflection,  both  of  which  could  modify  the  geometry  locally  (curved 
wall  and  surface),  but  are  unlikely  to  change  the  basic  structure  and  scaling. 

The  structure  of  the  thermal  and  flow  fields  will  be  examined  both  through  a  detailed 
scaling  analysis  to  determine  the  dependence  on  the  parameters,  and  through  numerical 
simulation  using  two  different  methods. 

2.  PROBLEM  STATMENT 

A  pool  of  incompressible  Newtonian  fluid  is  bounded  on  the  left  by  a  vertical  solid 
wall,  with  the  upper  portion  of  the  wall  (to  depth  d)  maintained  at  a  cold  temperature  Tc. 
while  the  rest  of  the  wall  is  at  the  hot  ambient  temperature  Th  of  undisturbed  fluid  far  from 
the  corner.  (See  figure  1.)  Above  the  horizontal  free  surface  of  the  liquid  is  an  inviscid, 
nonconducting  gas.  Surface  tension  is  assumed  strong  enough  to  keep  the  free  surface  flat 
(small  Capillary  number),  but  with  surface  tension  variations  due  to  a  linear  dependence 
on  temperature.  The  resulting  flow  is  assumed  to  be  two-dimensional  and  steady. 

Then  the  equations  governing  the  thermocapillary  convection  in  the  cold  corner  are 
conservation  of  mass,  momentum,  and  energy: 


V    u  =  0 

(2.1) 

pu 

•  Vu  =  -Vp  +  // V2u 

(2.2) 

pcpu 

•  V  T  =  k  V2  T 

(2.3) 

with  the  boundary  conditions: 

at     y  =  0  : 

Ty      =      0,                  V      =      0,                  flUy 

=  lTx 

(2.4a,b,c) 

at     x  —  0  : 

T  =  f  Tc,       y<d 
\Th,        y>d' 

U  =  V  = 

=  0 

(2.5a,b,c) 

as     x.  y  — >  oo  : 

T  — »  T/,.,       u,  v  — »  0 

(2.6a.b,c) 

Here  u  is  the  velocity  vector  with  components  u  and  v  in  the  x  (horizontally  rightward) 
and  y  (vertically  downward)  directions,  p  is  pressure,  T  is  temperature,  p  is  density,  p,  is 
viscosity,  cp  is  specific  heat,  k  is  thermal  conductivity,  and  7  (assumed  constant  and  neg- 
ative) is  the  derivative  of  the  surface  tension  with  respect  to  temperature.  The  boundary 
conditions  specify  that  the  wall  is  piecewise  isothermal  with  no  fluid  slip,  and  the  flat  free 
surface  is  thermally  insulated,  with  thermocapillary  forcing. 

The  equations  can  be  nondimensionalized  by  scaling  lengths  by  d.  temperature  dif- 
ferences by  AT  =  Th  —  Tc.  and  velocities  by  u3  =  7  AT/ p.  The  resulting  dimensionless 
equations  are: 

V-u  =  0  (2.7) 

RuVu=  -Vp  +  V2u  (2.8) 

Mu-VT  =  V2T  (2.9) 


i?u 

•  V 

U) 

= 

V   a; 
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= 

-V2 

* 

u  = 

, 

v  =  ■ 

-*« 

with  the  boundary  conditions: 

at     y  =  0  :      Ty  =  0.     u  =  0,     i^  =  Tx  (2.10a,b,c) 

f  -1.         v  <  1 
at     .t  =  0  :      T  =  <^  .      u  =  v  =  0  (2.11a,b,c) 

10,  y  >  1 

as     .t.  y  — ►  oo  :      T  — >  0.      k,  v  — >  0  (2.12a,b,c) 

where  u,  etc.  from  here  on  denote  the  dimensionless  quantities.  The  two  dimcnsionless 
parameters  are  the  Marangoni  number  M  =  usdj k  and  the  Reynolds  number  R  =  usd/v. 
where  k  =  k/ pcp  is  the  thermal  difFusivity  and  v  =  \if p  is  the  kinematic  viscosity.  Their 
ratio  gives  the  Prandtl  number:  P  =  u/k  =  M/R. 

For  the  numerical  solutions  below,  it  is  convenient  to  eliminate  the  pressure  by  adopt- 
ing a  stream-function/voiticity  formulation  for  the  flow: 

(2.13) 

(2.14) 

(2.15a,b) 

where  $  is  the  stream  function  and  to  is  the  vorticity,  with  the  flow  boundary  conditions 

at     y  =  0  :      *  =  *x  =  0,     u  =  -Tx  (2.16a,b,c) 

at     x  =  0  :      *  =  *x  =  Vy  -  0  (2.17a,b,c) 

as     x.y ->  oc  :      *.  u  -*  0  (2.18a.b) 

3.  SCALING  ANALYSIS  AND  REGIMES  OF  BEHAVIOR 

The  structure  of  the  thermal  and  flow  fields  can  take  on  different  forms,  depending  on 
the  values  of  the  two  governing  parameters,  the  Marangoni  number  M.  which  measures  the 
importance  of  thermal  convection  relative  to  thermal  diffusion,  and  the  Prandtl  number 
P,  which  is  the  ratio  of  viscous  to  thermal  diffusion,  a  material  property.  (Equivalently, 
one  could  use  Reynolds  number  R  =  M/P  as  the  second  parameter,  the  ratio  of  inertial  to 
viscous  forces.)  Here  we  derive  the  appropriate  dimensionless  scales  for  the  four  different 
asymptotic  regimes  of  behavior. 

To  examine  the  dominant  balances  in  the  cold  corner,  three  scales  suffice:  the  hori- 
zontal length  scale  I  for  the  thermal  gradient  along  the  surface,  the  vertical  viscous  length 
scale  8  for  the  velocity  shear  at  the  surface  (which  turns  out  to  be  the  same  scale  as  for 
the  velocity  shear  on  the  wall),  and  the  velocity  scale  U  for  flow  along  the  surface.  The 
vertical  thermal  length  scale  is  determined  by  the  boundary  condition  on  the  wall,  and  so 
is  0(1).  The  thermocapillary  stress  condition  (2.10c)  scales  as: 


so  that  U  ~  6/1.    In  the  energy  and  vorticity  equations  (2.9  and  2.13),  the  terms  for 
convection  in  each  direction  scale  the  same,  but  not  so  for  diffusion,  and  by  (2.1G),  the 
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1 

surface  vorticity  scale  is  1//,  so  (2.9)  and  (2.13)  scale  as: 

(3.2) 
(3-3) 

For  small  enough  M .  thermal  convection  is  negligible,  implying  /  ~  1,  and  the  thermal 
field  is  essentially  conductive,  decoupled  from  the  flow.  But  for  large  enough  M  convection 
becomes  important,  and  the  strong  surface  flow  toward  the  wall  compresses  the  thermal 
gradient  along  the  surface,  which  in  turn  strengthens  the  local  driving  force  for  the  flow. 
This  reduces  the  horizontal  thermal  length  scale  /  to  the  point  that  thermal  diffusion  away 
from  the  wall  balances  convection  toward  the  wall,  so  the  local  effective  Marangoni  number 
is  order  unity:  Meff  =  MUl  ~  1.  Then  the  externally  imposed  length  scale  (dimensional 
d  above)  is  no  longer  directly  relevant  to  the  compressed  cold  corner  region.  (In  this  case, 
the  local  importance  of  inertia  is  better  indicated  by  whether  viscous  or  thermal  diffusion 
is  more  efficient,  i.e..  by  P  rather  than  R.) 

Similarly,  for  small  enough  R.  inertia  is  negligible  everywhere,  implying  8  ~  /,  and  the 
flow  is  dominated  by  viscous  forces.  For  large  enough  R  incrtial  forces  become  dominant 
and  viscous  effects  are  confined  to  boundary  layers  of  thickness  8  <C  /  along  the  surface  and 
the  wall,  where  the  local  effective  Reynolds  number  Reff  =  RU82  /I  ~  1.  (Both  layers  are 
of  comparable  thickness  because  the  pressure  field  outside  the  layers  has  the  same  length 
scale  in  both  directions.)  The  above  gives  the  scaling  for  each  regime. 

When  the  thermal  field  is  conductive  and  the  flow  field  dominated  by  viscous  forces 
(M  <C  1  and  R  <  1.  or  P  >  M),  all  three  scales  are  of  order  unity:  I  ~  1,  8  ~  1,  U  ~  1. 
Thus  in  this  case  (only),  the  scaling  used  in  the  nondimensionalization  is  appropriate 
everywhere.  Within  this  regime,  the  solution  is  fully  two-dimensional  with  no  fine  structure 
and  is  nearly  independent  of  the  parameters. 

For  the  conductive  case  with  inertial  flow,  the  additional  resistance  of  inertia  reduces 
both  the  velocity  scale  and  the  viscous  length  scale:  U  ~  8  ~  R~s  (while  /  ~  1  still). 
This  reduced  velocity  also  reduces  the  effective  Marangoni  number,  such  that  this  regime 
applies  when  M  <C  R K  or  M  <  P~K  with  R  >  1,  or  P  <  M.  (Note  that  this  gives 
the  same  boundary  layer  scaling  as  Zcbib  et  al.,  1985,  except  they  made  an  error  on  the 
passive  wall  layers,  as  pointed  out  by  Chen,  1987.)  Here  the  vorticity  generated  by  the 
shear  stresses  on  the  surface  and  the  wall  are  confined  to  the  thin  boundary  layers. 

When  thermal  convection  is  important  but  inertia  is  not  (M  ^>  1  and  P  ^>  1),  surface 
thermal  varitions  are  compressed  to  a  narrow  region,  beyond  which  the  thermocapillary 
forcing  is  small,  so  8  ~  /  ~  M_1  and  U  ~  1.  However,  the  strong  inward  flow  along 
the  surface  turns  downward  and  away  from  the  no-slip  wall  (and  weakens  rapidly  with 
distance),  such  that  no  thermal  boundary  layer  is  formed  on  the  wall:  rather,  vertical  and 
horizontal  variations  are  comparable. 

The  most  important  regime  for  materials  processing  is  where  thermal  convection  is 
important  and  P<1,  the  latter  being  generally  true  for  metals.  In  this  case,  within  the 
compressed  thermal  region  there  are  thin  viscous  boundary  layers  along  the  surface  and 
wall.  Then  /  ~  M-1P~2 ,  8  ~  M-1,  and  U  ~  ?5 ,  i.e.,  the  additional  resistance  of  inertia 


decreases  the  velocity  scale  and  thus  increases  the  thermal  scale  by  a  factor  of  vP  relative 
to  the  purely  viscous  case.  Again  the  reduced  velocity  changes  the  thermal  convection 
scaling,  and  large  Marangoni  number  here  means  M  ^>  P~* . 

The  approximate  divisions  between  the  four  asymptotic  regimes  are  shown  in  figure 
2.  For  a  material  of  small  P,  as  M  is  increased  from  zero,  at  first  the  temperature  field  is 
conductive  and  the  flow  dominated  by  viscous  forces,  then  the  flow  becomes  primarily  iner- 
tial  and  viscous  boundaiy  layers  form,  and  finally  thermal  convection  becomes  dominant, 
shrinking  all  local  length  scales  in  the  corner. 

4.  VISCOUS  CORNER  REGION 

There  is  a  region  in  the  corner,  for  any  M  and  P.  where  viscous  stresses  from  the  wall 
limit  the  flow  and  both  inertia  and  thermal  convection  arc  negligible,  so  the  temperature 
is  a  linear  function  along  the  surface.  Locally  the  thermocapillary  stress  is  constant,  and 
the  flow  is  given  by  a  similarity  solution  (Moffatt,  19G4,  although  the  published  version  is 
incorrect). 

If  the  flat  free  surface  makes  an  angle  a  with  the  solid  wall,  then  a  constant  unit 
surface  shear  stress  toward  the  corner  gives 


#(r,0)  = 


r2  [(sin  2a  -  2a)(cos20  -  1)  -  (cos  2a  -  l)(sin20  -  20)] 
4  sin  2a  —  2a  cos  2a 


(4.1) 


where  r  and  0  are  polar  coordinates,  with  0  increasing  from  the  wall  to  the  free  surface 
(see  figure  1).  In  the  special  case  here,  a 


j,  and 


$(r,0)  =  r: 


u  =  r 


-(1  -cos 20)4-  —(sin 20 -20) 
4  2tt 


-sin  20  +  -(cos  20-  1) 
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where  f  and  0  are  unit  vectors  in  the  coordinate  directions,  and  po  is  some  reference 
pressure.  Figure  3  shows  the  streamlines  for  this  similarity  flow.  The  zero-vorticity  contour 
extends  from  the  corner  at  a  45°  angle,  dividing  the  negative  vorticity  on  the  surface  from 
the  positive  vorticity  on  the  wall. 

This  is  the  form  of  the  flow  in  the  cold  corner  on  the  smallest  length  scale,  where 
all  the  above  flow  quantities  would  be  multiplied  by  Tx(0,0)  (which  scales  as  1//).  The 
velocity  grows  linearly  with  distance  r  from  the  corner  and  the  local  length  scale  is  r,  which 
can  be  used  to  estimate  the  range  over  which  the  similarity  solution  is  applicable.  For  the 
two  viscous-flow  regimes  mentioned  above,  the  linear-temperature  approximation  requires 
r  <C  I,  so  for  the  conductive  regime  r  <C  1  while  for  the  convective  regime  r  <C  M-1. 


For  the  two  inertial-flow  regimes,  the  more  restrictive  condition  is  that  locally  inertia  is 
negligible,  and  the  local  velocity  scale  is  r/7,  so  we  must  have  r2  <C  l/R;  in  the  conductive 
regime  this  implies  r  <C  i?".  while  in  the  convectivc  regime  the  range  is  r  <C  P1  M_1,  or 
r/KPt 

This  gives  an  estimate  of  the  resolution  required  for  a  numerical  model  to  resolve  all 
the  details  of  the  cold  corner  flow.  Because  of  the  corner  singularity  in  u  and  p,  spectral 
methods  would  be  inappropriate;  instead,  finite-difference  or  finite-element  methods  could 
be  used.  Then  if  the  first  grid  point  is  in  the  similarity  range,  no  details  of  the  velocity  field 
in  the  corner  will  be  lost.  In  addition,  the  similarity  form  may  be  useful  as  a  "matching" 
type  boundary  condition  for  the  singularity  at  the  origin. 

5.  GREEN'S  FUNCTION  METHOD  FOR  VISCOUS  CASE:  R  ->  0 

When  the  Reynolds  number  is  sufficiently  small,  inertia  is  negligible  throughout  the 
flow  field,  and  to  lowest  approximation  the  flow  is  governed  by  the  biharmonic  equation: 

V2(V2*)  =  0  (5.1) 

Then  the  flow  everywhere  depends  only  on  the  instantaneous  thermal  gradient  along  the 
surface  (even  if  the  flow  is  unsteady).  This  allows  the  flow  field  to  be  represented  using 
the  Green's  function  for  a  point  force  near  a  rigid  wall,  directed  toward  the  wall  (see,  e.g., 
Hasimoto  &  Sano.  1980,  or  Blake,  1971). 

The  G™een1s  function  for  the  stream  function  due  to  a  unit  (dimensionless)  force  in 
the  negative  x  direction  applied  on  the  surface  at  the  point  (£,  0)  is: 


Z7T 


V  I" 


(5.2a) 


where 

r+  =  ^(x  +  02+y2-      r_  =  ^(x-tf+y2  (5.2b) 

The  flow  field  due  to  such  a  point  force  at  £  =  1  is  shown  in  figure  4.  (Note  that 
all  distances  scale  with  £.)  The  far  field  (r  ^>  £)  can  be  used  as  an  approximation  to  the 
viscous  thermocapillary  flow  at  large  enough  distances  that  the  distributed  forcing  of  the 
thermal  gradient  can  be  replaced  by  a  point  force  near  the  wall.  It  can  be  shown  that  in 
the  far-field  the  zero  vorticity  contour  approaches  an  angle  of  j  from  the  surface.  As  r 
increases,  the  velocity  decays  rapidly  as  u  ~  0(r~2). 

This  rapid  decay  of  velocity  has  the  consequence  that  no  thermal  boundary  layer 
forms  on  the  wall.  While  this  result  only  applies  directly  when  inertia  is  negligible,  one 
would  expect  that  the  same  order  of  decay  applies  to  the  inertial  case  in  the  far  field,  by 
analogy  with  the  Squire-Landau  jet  (Landau,  1944,  and  Squire,  1951)  where  the  velocity 
due  to  an  isolated  point  force  decays  like  u  ~  0(r_1)  for  all  Reynolds  numbers.  Note  also 
that  near  the  wall  (except  near  the  surface),  the  flow  is  away  from  the  wall.  Hence,  though 
the  surface  temperature  gradient  may  be  highly  compressed  from  vigorous  convection, 
there  is  no  thermal  boundary  layer  along  the  wall,  i.e.,  vertical  and  horizontal  temperature 
variations  are  expected  to  be  of  the  same  order.  This  is  in  marked  contrast  with  the 
seemingly  comparable  case  of  a  hot  corner.   There  the  flow  tends  to  expand  the  surface 


thermal  gradient  over  long  distance,  so  the  forcing  is  similar  to  a  distributed  line  force, 
and  thermal  boundary  layers  form  both  on  the  wall  and  the  surface  (Cowley  and  Davis. 
1983).  This  difference  between  concentrated  forcing  and  distributed  forcing  changes  the 
fundamental  character  of  the  flow. 

Using  the  Green's  function,  then  the  stream  function,  velocity  components,  and  vor- 
ticity  in  the  viscous  thermocapillary  flow  are  given  by: 
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using  the  shortened  notation  T'(x)  =  Tx(x,0). 

To  calculate  the  steady  flow  for  various  values  of  Marangoni  number  (with  i?  =  0).  the 
time-dependent  equation  for  the  temperature  field  was  explicitly  integrated  in  time  until 
the  steady  state  was  reached.  A  finite  square  domain  on  a  uniform  grid  was  used,  where  on 
the  artificial  boundaries,  diffusion  across  the  boundaries  is  neglected  and  where  convection 
is  inward  the  fluid  outside  is  assumed  to  be  at  zero  temperature.  The  effects  of  these 
artificial  boundaries  on  the  cold  corner  are  presumably  small  if  the  boundaries  are  several 
units  away,  since  the  thermal  field  decays  quickly  with  distance.  Central  finite  differences 
were  used,  with  upwind  differencing  for  the  convective  terms.  The  velocity  at  each  point 
was  evaluated  from  the  above  integrals  (5.3b.c).  using  an  analytic  approximation  around 
the  singularity  at  (x,  0).  and  elsewhere  using  the  trapezoidal  rule  with  first-order  differences 
for  the  thermal  gradient,  as  central  differences  here  resulted  in  a  numerical  instability. 

The  results  for  a  range  of  M  are  shown  in  figure  5.  The  transition  from  the  primarily 
conductive  regime  to  the  convective  regime  as  M  increases  is  apparent.  As  the  surface 
gradient  becomes  more  compact,  the  strong  downward  convection  away  from  the  corner 
extends  the  thermal  region  into  a  sort  of  broad  diagonal  plume;  no  thermal  boundary  layer 
forms  even  for  large  M.  (Where  the  isotherms  intersect  the  artificial  boundaries,  slight 
local  distortion  due  to  the  artificial  boundary  conditions  is  apparent.)  For  large  M,  the 
flow  outside  the  cold  corner  is  qualitatively  similar  to  the  flow  from  a  point  source  (see 
figure  4). 

G.  NUMERICAL  SOLUTIONS  FOR,  FULL  SYSTEM 

When  inertia  is  not  negligible,  the  full  coupled  nonlinear  equations  must  be  solved. 
Unlike  the  Green's  function  approach,  this  requires  flow  boundary  conditions  at  the  arti- 
ficial boundaries  of  the  computational  domain;  here  these  boundaries  are  assumed  to  be 
impermeable  and  shear-free,  as  well  as  isothermal.  These  conditions  constrain  both  the 
thermal  and  flow  fields  (compared  to  the  Green's  function  method),  enhancing  recircula- 
tion and  preventing  long  thermal  plumes.  Still,  with  the  artificial  boundaries  several  units 
away,  their  effects  on  the  cold  corner  are  expected  to  be  small. 


Again,  the  numerical  method  involves  stepping  the  unsteady  equations  forward  in  time 
until  steady  state  is  reached.  At  each  time  step,  the  convection-diffusion  equations  for  tem- 
perature (2.9)  and  vorticity  (2.13)  are  solved  by  the  Alternating-Direction  Implicit  (ADI) 
scheme,  where  the  convective  terms  are  evaluated  by  the  Eulerian-Lagrangian  Method 
(ELM.  see  Cheng  et  al..  1984)  using  the  velocity  field  from  the  old  stream  function  and 
upwind  bilinear  interpolation.  (The  ADI  method  avoids  diffusive  numerical  instability,  so 
the  time  step  is  only  limited  by  convective  stability.)  After  several  time  steps,  the  Poisson 
equation  for  the  stream  function  (2.14)  is  solved  by  Gauss-Seidel  iteration  with  Succes- 
sive Over- Relaxation  (SOR).  Steady  state  is  assumed  when  the  pointwise  RMS  change  in 
stream  function  is  below  a  certain  tolerance. 

A  non-uniform  Cartesian  grid  was  employed.  (The  program  allows  arbitrary  spacing 
of  points  in  each  direction.)  Using  a  three-point  difference  scheme,  only  first-order  accuracy 
is  possible  for  the  second  derivatives;  the  differencing  employed  becomes  second-order  in 
the  limit  of  uniform  spacing.  Specifically,  the  following  one-dimensional  difference  formulae 
were  used  (derived  from  Taylor  series): 

r,  =  ~<fc+  ,      |    dx+  -  dx.  dx_ 

dx-(dx+  +  dx-)  dx+dx-  dx+{dx+  +  dx-) 

+  0{f'"dx+dx_)  (6.1a) 

2  -2  2 

f"  — r    _i t  _|_ r 

dx-(dx+  +  dx.)  dx+dx.  dxjr{dx+  +  dx—) 

+  0{f"'{dx+-dx-))  (6.1b) 

where  dx.  and  dx+  are  the  distances  to  the  grid  points  below  and  above  the  current  point, 
with  /_  and  /+  the  corresponding  values  of  the  function. 

The  grid  spacing  in  each  direction  was  chosen  to  have  a  closely  packed  region  of 
uniform  spacing  by  the  surface  or  wall,  a  widely  spaced  region  of  uniform  spacing  out  near 
the  artificial  boundaries,  and  in  between  a  region  of  smoothly  (exponentially)  changing 
spacing.  This  was  generated  by  applying  the  following  function  to  a  uniform  grid  (in  £. 
say): 

lnr£/A 

\  (6-6)r 

where  r  is  a  parameter  giving  the  ratio  of  outer  to  inner  spacing,  £i  and  £2  are  parameters 
delineating  the  three  regions.  D  =  (r  —  1)(£2  —  £1)  +  h\r(r  +  £1  —  ^£2),  and  xmax  is  the 
position  of  the  artificial  boundary. 

It  was  found  that,  even  when  the  time  step  easily  satisfied  convective  stability  require- 
ments, nonetheless  instabilities  developed  in  the  vorticity  near  the  wall.  Several  different 
formulations  for  the  wall  boundary  condition  on  vorticity  were  tried,  to  no  avail.  However, 
by  under-relaxing  the  changes  in  wall  vorticity  only  during  the  initial  adjustment  period, 
the  instability  was  eliminated.  The  formula  used  to  calculate  the  vorticity  at  the  wall 
(without  under-relaxation)  from  the  stream  function  is 

u0  =  2(*2:  '  *2)      +0(*'"x2)  (6.3) 
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where  subscripts  0,  1,  and  2  refer  to  the  wall  and  the  first  two  grid  points.  While  only 
first-order,  this  formulation  is  independent  of  \£o  and  thus  avoids  re-using  the  boundary 
condition  for  \I>. 

Results  are  shown  in  figures  6  and  7  for  P  =  1  and  P  =  0.01,  for  a  variety  of  M .  The 
conductive,  viscous  regime  is  represented  by  the  case  M  =  0.01,  P  =  1  (figure  6a).  For 
P  =  1.  with  increasing  M,  the  surface  thermal  gradient  becomes  compressed  (similar  to 
the  P  =  oc  case  computed  by  Green's  functions)  along  with  the  vorticity  on  the  surface. 
However,  inertia  is  no  longer  negligible,  and  so  the  flow  down  the  wall  has  less  of  an 
outward  component.  Also,  the  artificial  boundaries  modify  the  ''plume,"  keeping  it  from 
the  boundaries  and  turning  it  upward  due  to  the  recirculation.  For  P  —  0.01  (as  typical  in 
liquid  metals),  with  increasing  M  inertial  effects  confine  the  vorticity,  forming  clear  viscous 
boundary  layers,  before  thermal  convection  becomes  strong.  Also,  a  counter-rotating  cell 
forms  in  the  lower  part  of  the  domain.  As  thermal  convection  compresses  the  surface 
gradient,  the  surface  viscous  boundary  layer  remains  limited  to  the  cold  corner.  Note  that 
for  large  M  and  small  P,  the  thermal  length  scale  is  small  compared  to  the  overall  domain, 
and  the  viscous  length  scale  is  small  compared  to  the  thermal  scale,  imposing  severe  (local) 
resolution  requirements  on  any  numerical  model. 

The  numerical  results  for  P  =  1  and  P  =  0.01  are  compared  with  the  scaling  analysis 
in  figure  8.  The  thermal  length  scale  (or  rather  1//)  is  estimated  from  the  thermal  gradient 
Tx  at  the  wall  (based  on  the  first  grid  point  along  the  surface).  The  velocity  scale  (U) 
is  estimated  from  the  maximum  velocity  umax  at  a  grid  point  on  the  surface  (though  the 
actual  maximum  might  be  expected  to  fall  between  grid  points).  The  viscous  length  scale 
(6)  is  approximated  by  the  position  xmax  of  the  grid  point  with  the  maximum  velocity. 
For  P  =  1,  the  transition  is  apparent  from  the  conductive  viscous  regime  (where  /,  U, 
and  6  are  nearly  independent  of  M)  to  the  convective  viscous  regime  (where  /  oc  M~l  and 
6  oc  M_1,  with  U  roughly  constant).  For  P  =  0.01,  two  transitions  are  seen,  from  the 
conductive  viscous  regime  to  the  conductive  inertia!  regime  (/  nearly  constant,  U  oc  M-1'3. 
and  8  oc  M_1//3).  to  the  convective  inertia!  regime  (I  oc  M_1,  6  oc  M~l .  U  nearly  constant). 

7.  DISCUSSION  AND  CONCLUSIONS 

The  practical  importance  of  thermocapillary  convection  in  materials  processing,  along 
with  the  complications  inherent  in  typical  processes  (e.g.  curved  interfaces,  phase  change, 
etc.),  insure  that  numerical  simulations  will  remain  one  of  the  main  theoretical  tools  for 
understanding  such  systems.  This  work  predicts,  a  priori,  the  resolution  requirements  for 
such  numerical  models  to  accurately  represent  the  high  heat  transfer  and  rapid  velocity 
variations  in  the  cold  corner  region. 

The  structure  of  the  corner  depends  on  two  dimensionless  parameters  indicating  the 
driving  force  for  convection  and  relative  importance  of  viscosity:  the  Marangoni  number 
M  (based  on  the  overall  temperature  difference,  overall  thermal  length  scale,  and  material 
properties)  and  the  Prandtl  number  P  (a  material  property).  Hence  there  are  four  asymp- 
totic regimes  (shown  in  figure  2)  depending  on  whether  thermal  convection  and  inertial 
forces  are  locally  important.  For  large- P  materials  (e.g.,  organics),  the  corner  flow  is  dom- 
inated by  viscous  forces,  and  for  large  M  the  surface  thermal  variations  are  compressed, 
i.e.,  the  local  length  scale  is  decreased  (and  heat  transfer  increased).  For  small-P  materi- 
als (e.g.,  metallics),  inertia  becomes  important  before  thermal  convection  with  increasing 
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M,  and  thin  viscous  boundary  layers  form  within  the  thermal  region.  When  M  is  large 
(compared  to  P-1'2),  three  levels  of  length  scales  must  be  resolved  (overall,  thermal,  and 
viscous),  a  severe  requirement  on  numerical  models. 

The  scaling  was  derived  from  a  simple  problem  designed  to  isolate  the  feedback  mech- 
anism of  the  cold  corner.  Numerical  calculations  by  two  methods  for  closely  related  prob- 
lems (where  the  corner  is  necessarily  less  isolated  due  to  the  finite  domain)  illustrate  the 
changing  structure  of  the  cold  corner  and  show  the  details  of  the  transitions  between  the 
asymptotic  regimes.  The  numerical  results  are  consistent  with  the  scaling  analysis. 

One  surprising  result,  in  contrast  with  the  hot  corner  problem  of  Cowley  and  Davis 
(1983),  is  that  no  thermal  boundary  layers  form.  This  difference  is  due  to  the  forcing  being 
limited  to  a  relatively  concentrated  region  in  the  cold  corner,  while  for  the  hot  corner  the 
forcing  is  distributed  over  an  extended  region  (the  horizontal  thermal  variations  being 
extended  by  convection). 

To  compare  with  the  numerical  results  of  Zebib  et  al.  (1985)  for  P  =  1,  note  the 
different  boundary  conditions:  their  domain  had  a  (flat)  free  surface,  a  hot  and  a  cold 
no-slip  wall,  and  an  insulated  no-slip  bottom.  Then  even  when  the  cold  corner  was  highly 
compressed,  there  were  still  significant  thermal  variations  along  most  of  the  surface  due  to 
the  hot  wall,  so  the  overall  scaling  is  similar  to  the  conductive  inertia]  case  here.  And  ap- 
parently even  their  cold  corner  was  modified  by  the  external  bulk  flow,  for  their  maximum 
vorticity  scaled  as  woe  M2'3,  whereas  here  u  <x  M.  (The  resolution  requirements  derived 
here  can  thus  be  considered  conservative.)  Hence  one  important  question  remaining  is 
under  what  conditions  can  the  cold  corner  be  considered  locally  determined  (as  it  is  here). 

In  real  materials  processes,  the  surface  is  free  to  deflect,  and  the  position  and  shape  of 
the  solid-liquid  interface  depends  on  the  thermal  field.  These  effects  greatly  complicate  the 
problem  geometrically,  yet  the  dominant  dynamic  balance  should  remain  that  considered 
here.  Future  work  will  investigate  these  effects  in  the  cold  corner. 
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FIGURE  CAPTIONS 

Figure  1.  Problem  Formulation:  a  liquid  quarter-space  is  bounded  above  by  a  flat  free 
surface  subject  to  thermocapillary  forcing,  and  is  bounded  on  the  left  by  a  rigid  vertical 
wall,  at  temperature  Tc  to  depth  d  and  at  the  warmer  temperature  Th  below,  which  is  also 
the  ambient  temperature  of  the  undisturbed  fluid  far  away. 

Figure  2.  Asymptotic  Regimes:  the  four  different  asymptotic  behaviors  correspond  to 
whether  the  thermal  field  is  controlled  primarily  by  convection  or  by  conduction,  and 
whether  the  flow  is  dominated  by  viscous  or  inertial  forces. 

Figure  3.  Viscous  Corner  Flow:  streamlines  (solid,  black)  and  vorticity  contours  (dashed, 
light  gray)  for  the  self- similar  solution  valid  where  the  surface  temperature  is  approximately 
linear  and  inertia  is  negligible. 

Figure  4.  Green's  Function:  streamlines  (solid,  black)  and  vorticity  contours  (dashed,  light 
gray)  for  the  non-inertial  flow  due  to  a  point  force  on  the  surface  at  (1,0),  directed  toward 
the  wall.  The  flow  is  recirculating  and  decays  rapidly  with  distance. 

Figure  5.  Numerical  Results  for  i?  — >  0  (P  — >  oc):  velocity  vectors  and  isotherms  (solid, 
dark  gray)  of  numerical  solutions  using  the  Green's  function  formulation,  for  a  range  of 
Marangoni  numbers.  On  the  artificial  boundaries  (bottom  and  right),  normal  thermal 
diffusion  is  neglected  and  incoming  fluid  is  assumed  isothermal,  (a)  M  =  10;  (b)  M  —  30: 
(c)  M  =  100;  (d)  M  =  300:  (e)  M  =  1000. 

Figure  6.  Numerical  Results  for  P  =  1:  isotherms  (solid,  dark  gray),  streamlines  (thin, 
black)  and  vorticity  contours  (dashed,  light  gray),  for  the  numerical  solution  of  the  full 
coupled  system.  The  artificial  boundaries  (bottom  and  right)  are  assumed  isothermal, 
impenetrable,  and  shear-free,  (a)  M  =  0.01:  (b)  M  =  10;  (c)  M  =  30;  (d)  M  =  100;  (e) 
M  =  300:  (f )  M  =  1000:  (g)  detail  of  M  =  1000. 

Figure  7.  Numerical  Results  for  P  =  0.01:  see  previous  caption,  (a)  M  =  1;  (b)  M  =  10; 
(c)  M  =  100;  (d)  detail  of  M  =  100:  (e)  M  =  1000;  (f)  detail  of  M  =  1000:  (g)  M  = 
10.000;  (h)  detail  of  M  =  10,000. 

Figure  8.  Summary  of  Numerical  Scales:  wall  temperature  gradient  Tx  (diamonds),  maxi- 
mum velocity  umax  (triangles),  and  position  xmax  of  maximum  velocity  (stars),  as  functions 
of  Marangoni  number  M,  from  numerical  solutions  of  full  system.  Lines  of  slope  1,-1.  and 
-1/3  are  shown  for  comparison  with  scaling  analysis,  (a)  P  =  1;  (b)  P  =  0.01. 
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